##### ##################################################################
#####
#####   Input: Clean experimental data
#####   Output: Analysis and figures for the supplementary materials
#####

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
rm(list=ls())

library(ggplot2)
library(data.table)
library(questionr)
library(tidyverse)
library(dplyr)
library(broom)
library(patchwork)
library(ggthemes)
library(stargazer)
library(knitr)

load("data/reshaped_data.Rdata")

scaleFUN <- function(x) sprintf("%.2f", x)

# Appendix -- Tables 2 & 3 ---------------------------------------------------------------------
## Unconditional trait preferences by crisis type
competence <- lm(Q2_1 ~ Condition, data = stacked_data, weight = GBw8)
honest <- lm(Q2_2 ~ Condition, data = stacked_data, weight = GBw8)
benevolent <- lm(Q2_3 ~ Condition, data = stacked_data, weight = GBw8)
authentic <- lm(Q2_4 ~ Condition, data = stacked_data, weight = GBw8)
caring <- lm(Q2_5 ~ Condition, data = stacked_data, weight = GBw8)
assertive <- lm(Q2_6 ~ Condition, data = stacked_data, weight = GBw8)
risk <- lm(Q2_7 ~ Condition, data = stacked_data, weight = GBw8)
collaborative <- lm(Q2_8 ~ Condition, data = stacked_data, weight = GBw8)
empathy <- lm(Q2_9 ~ Condition, data = stacked_data, weight = GBw8)

stargazer(
  competence, honest, benevolent, authentic, caring, 
  font.size = "scriptsize",
  dep.var.labels = c("Competent", "Honest", "Benevolent", "Authentic", "Caring"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Economic Crises", "Health Crises", "Security Crises"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Trait Preferences by Crisis Type")

stargazer(
  assertive, risk, collaborative, empathy,
  font.size = "scriptsize",
  dep.var.labels = c("Assertive", "Risk-Taker", "Collaborative", "Empathetic"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Economic Crises", "Health Crises", "Security Crises"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Traits by Sex of Respondent - All Crisis Types", label="tab:hyp1-3_women")

# Appendix -- Tables 4 & 5 ---------------------------------------------------------------------
## Unconditional trait preferences by respondent gender
competence <- lm(Q2_1 ~ gender, data = stacked_data, weight = GBw8)
honest <- lm(Q2_2 ~ gender, data = stacked_data, weight = GBw8)
benevolent <- lm(Q2_3 ~ gender, data = stacked_data, weight = GBw8)
authentic <- lm(Q2_4 ~ gender, data = stacked_data, weight = GBw8)
caring <- lm(Q2_5 ~ gender, data = stacked_data, weight = GBw8)
assertive <- lm(Q2_6 ~ gender, data = stacked_data, weight = GBw8)
risk <- lm(Q2_7 ~ gender, data = stacked_data, weight = GBw8)
collaborative <- lm(Q2_8 ~ gender, data = stacked_data, weight = GBw8)
empathy <- lm(Q2_9 ~ gender, data = stacked_data, weight = GBw8)

stargazer(
  competence, honest, benevolent, authentic, caring, 
  font.size = "scriptsize",
  dep.var.labels = c("Competent", "Honest", "Benevolent", "Authentic", "Caring"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Woman Respondent"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Traits by Sex of Respondent - All Crisis Types", label="tab:hyp1-3_women")

stargazer(
  assertive, risk, collaborative, empathy,
  font.size = "scriptsize",
  dep.var.labels = c("Assertive", "Risk-Taker", "Collaborative", "Empathetic"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Woman Respondent"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Traits by Sex of Respondent - All Crisis Types", label="tab:hyp1-3_women")

## Appendix -- Tables 6 & 7 ---------------------------------------------------------------------
## Conditional trait preferences by crisis type and respondent gender
competence <- lm(Q2_1 ~ Condition*gender, data = stacked_data, weight = GBw8)
honest <- lm(Q2_2 ~ Condition*gender, data = stacked_data, weight = GBw8)
benevolent <- lm(Q2_3 ~ Condition*gender, data = stacked_data, weight = GBw8)
authentic <- lm(Q2_4 ~ Condition*gender, data = stacked_data, weight = GBw8)
caring <- lm(Q2_5 ~ Condition*gender, data = stacked_data, weight = GBw8)
assertive <- lm(Q2_6 ~ Condition*gender, data = stacked_data, weight = GBw8)
risk <- lm(Q2_7 ~ Condition*gender, data = stacked_data, weight = GBw8)
collaborative <- lm(Q2_8 ~ Condition*gender, data = stacked_data, weight = GBw8)
empathetic <- lm(Q2_9 ~ Condition*gender, data = stacked_data, weight = GBw8)

stargazer(
  competence, honest, benevolent, authentic, caring, 
  font.size = "scriptsize",
  dep.var.labels = c("Competent", "Honest", "Benevolent", "Authentic", "Caring"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Economic","Health", "Security", 
                        "Woman", "Economic*Woman","Health*Woman", "Security*Woman"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Traits by Sex of Respondent - By Crisis Types", label="tab:hyp1-3_women")

stargazer(
assertive, risk, collaborative, empathetic,
  font.size = "scriptsize",
  dep.var.labels = c("Assertive", "Risk-Taker", "Collaborative", "Empathetic"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Economic","Health", "Security", 
                        "Woman", "Economic*Woman","Health*Woman", "Security*Woman"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Traits by Sex of Respondent - By Crisis Types", label="tab:hyp1-3_women")

# Appendix -- Table 8 ---------------------------------------------------------------------
## Marginal effect of gender (Women − Men) within each condition
marginal_effects_gender <- function(model, dv_label) {
  b <- coef(model); V <- vcov(model)
  conds <- c("Control","Economic","Health","Security")
  rows <- lapply(conds, function(cn) {
    if (cn == "Control") {
      est <- b["genderWoman"]
      se  <- sqrt(V["genderWoman","genderWoman"])
    } else {
      int <- paste0("Condition", cn, ":genderWoman")
      est <- b["genderWoman"] + b[int]
      se  <- sqrt(V["genderWoman","genderWoman"] +
                  V[int,int] + 2*V["genderWoman", int])
    }
    data.frame(
      Trait = dv_label,
      Condition = cn,
      MarginalEffect = est,
      SE = se,
      LowerCI = est - 1.96*se,
      UpperCI = est + 1.96*se,
      z = est / se,
      p = 2*pnorm(-abs(est / se))
    )
  })
  do.call(rbind, rows)
}

## Run on all models 
ME <- do.call(rbind, list(
  marginal_effects_gender(competence,   "Competence"),
  marginal_effects_gender(honest,       "Honest"),
  marginal_effects_gender(benevolent,   "Benevolent"),
  marginal_effects_gender(authentic,    "Authentic"),
  marginal_effects_gender(caring,       "Caring"),
  marginal_effects_gender(assertive,    "Assertive"),
  marginal_effects_gender(risk,         "Risk-taking"),
  marginal_effects_gender(collaborative,"Collaborative"),
  marginal_effects_gender(empathetic,   "Empathy")
))

ME_out <- within(ME, {
  MarginalEffect <- round(MarginalEffect, 2)
  SE             <- round(SE, 2)
  LowerCI        <- round(LowerCI, 2)
  UpperCI        <- round(UpperCI, 2)
  p              <- sprintf("%.2f", p)
})

kable(
  ME_out[, c("Trait", "Condition", "MarginalEffect", "SE", "LowerCI", "UpperCI")],
  format = "latex",
  booktabs = TRUE,
  digits = 3,
  row.names = FALSE,
  caption = "Marginal effect of gender (Women -- Men) on selection of traits by crisis type."
)

# Appendix -- Table 9 ---------------------------------------------------------------------
## Unconditional gendered traits by crisis type 
feminine <- lm(feminine_traits ~ Condition, data = stacked_data, 
   weight = GBw8)
masculine <- lm(masculine_traits ~ Condition, data = stacked_data, 
   weight = GBw8)
neutral <- lm(neutral_traits ~ Condition, data = stacked_data, 
   weight = GBw8)

stargazer(
  feminine, masculine, neutral,
  font.size = "scriptsize",
  dep.var.labels = c("Feminine traits", "Masculine traits", "Neutral traits"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
  covariate.labels = c("Intercept", "Economic","Health", "Security"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Preferable Gendered Traits - By Crisis Types", label="tab:hyp1-3_women")

# Appendix -- Table 10 ---------------------------------------------------------------------
## Conditional gendered traits by crisis type and respondent gender 
feminine <- lm(feminine_traits ~ Condition*gender, data = stacked_data, weight = GBw8)
masculine <- lm(masculine_traits ~ Condition*gender, data = stacked_data, weight = GBw8)
neutral <- lm(neutral_traits ~ Condition*gender, data = stacked_data, weight = GBw8)

stargazer(
  feminine, masculine, neutral,
  font.size = "scriptsize",
  dep.var.labels = c("Feminine traits", "Masculine traits", "Neutral traits"), 
  intercept.bottom = FALSE, 
  no.space=TRUE, keep.stat=c("n","rsq","adj.rsq"),
 covariate.labels = c("Intercept", "Economic","Health", "Security", 
                       "Woman", "Economic*Woman","Health*Woman", "Security*Woman"),
  header = FALSE, star.cutoffs = c(0.05, 0.01, 0.001),
  title="Preferable Gendered Traits - By Crisis Types and Respondent Sex", label="tab:hyp1-3_women")

# Appendix -- Table 11 ---------------------------------------------------------------------
## Marginal effect of gender (Women − Men) within each condition

marginal_effects_gender <- function(model, dv_label) {
  coefs <- coef(model)
  vcov_mat <- vcov(model)

  conditions <- c("Control", "Economic", "Health", "Security")
  results <- list()
  
  for (cond in conditions) {
    if (cond == "Control") {
      est <- coefs["genderWoman"]
      se <- sqrt(vcov_mat["genderWoman", "genderWoman"])
    } else {
      term_woman <- "genderWoman"
      term_int <- paste0("Condition", cond, ":genderWoman")
      est <- coefs[term_woman] + coefs[term_int]
      
      se <- sqrt(
        vcov_mat[term_woman, term_woman] +
        vcov_mat[term_int, term_int] +
        2 * vcov_mat[term_woman, term_int]
      )
    }
    
    lower <- est - 1.96 * se
    upper <- est + 1.96 * se
    
    results[[cond]] <- data.frame(
      Trait_Category = dv_label,
      Condition = cond,
      MarginalEffect = est,
      SE = se,
      LowerCI = lower,
      UpperCI = upper
    )
  }
  
  do.call(rbind, results)
}

## Run for all three models
me_feminine <- marginal_effects_gender(feminine, "Feminine traits")
me_masculine <- marginal_effects_gender(masculine, "Masculine traits")
me_neutral <- marginal_effects_gender(neutral, "Neutral traits")

## Combine into one table
me_all <- rbind(me_feminine, me_masculine, me_neutral)

## Round 
me_all$MarginalEffect <- round(me_all$MarginalEffect, 2)
me_all$SE <- round(me_all$SE, 2)
me_all$LowerCI <- round(me_all$LowerCI, 2)
me_all$UpperCI <- round(me_all$UpperCI, 2)



